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Abstract 

We present a practical algorithm for partially relaxing multiwell energy densities such as 
pertain to materials undergoing martensitic phase transitions. The algorithm is based on se- 
quential lamination, but the evolution of the microstructure during a deformation process is 
required to satisfy a continuity constraint, in the sense that the new microstructure should be 
reachable from the preceding one by a combination of branching and pruning operations. All 
microstructures generated by the algorithm are in static and configurational equilibrium. Ow- 
ing to the continuity constrained imposed upon the microstructural evolution, the predicted 
material behavior may be path-dependent and exhibit hysteresis. In cases in which there is a 
strict separation of micro and macrostructural lengthscales, the proposed relaxation algorithm 
may effectively be integrated into macroscopic finite-element calculations at the subgrid level. 
We demonstrate this aspect of the algorithm by means of a numerical example concerned with 
the indentation of an Cu-Al-Ni shape memory alloy by a spherical indenter. 
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1 Introduction 



Materials often are capable of adopting a multiplicity of crystal structures, or phases, the relative 
stability of which depends on temperature, the state of stress, and other factors. Under conditions 
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such that several phases are energetically favorable, e. g., at the transition temperature in marten- 
sitic materials, materials are often found to develop microstructure in nature or in the laboratory. 
A central problem in mechanics concerns the prediction of these microstructures and their effect 
on the effective or macroscopic behavior of materials, including such scaling properties and size 
effects as may result from their formation and evolution. When martensitic materials are modelled 
within the confines of nonlinear elasticity, the coexistence of phases confers their strain-energy 
density function a multiwell structure [TJ], ^5]] . The corresponding boundary value prob- 

lems are characterized by energy functions which lack weak sequential lower-semicontinuity, and 
the energy-minimizing deformation fields tend to develop fine microstructure [|7], ^3], [37p. 

There remains a need at present for efficient numerical methods for solving macroscopic 
boundary-value problems while simultaneously accounting for microstructure development at the 
microscale. One numerical strategy consists of attempting a direct minimization of a suitably 
discretized energy function. For instance, Tadmor et al [ |48| ] have applied this approach to the 
simulation of nanoindentation in silicon. The energy density is derived from the Stillinger- Weber 
potential by recourse to the Cauchy-Born approximation, and accounts for five phases of silicon. 
The energy functional is discretized by an application of the finite-element method. Tadmor et al 
[[49]] pioneering calculations predict the formation of complex phase arrangements under the inden- 
ter, and such experimentally observed features as an insulator-to-conductor transition at a certain 
critical depth of indentation. 

Despite these successes, direct energy minimization is not without shortcomings. Thus, anal- 
ysis has shown (see, e. g., [ |T7| ] for a review) that the microstructures which most effectively relax 
the energy may be exceedingly intricate and, consequently, unlikely to be adequately resolved by 
a fixed numerical grid. As a result, the computed microstructure is often coarse and biased by 
the computational mesh, which inhibits — or entirely suppresses — the development of many of the 
competing microstructures. By virtue of these constraints, the numerical solution is often caught 
up in a metastable local minimum which may not accurately reflect the actual energetics and de- 
formation characteristics of the material. 

In applications where there is a clear separation of micro and macro structural lengthscales, an 



alternative numerical strategy is to use a suitably relaxed energy density in calculations [ j20| , 
[10|, ^4], [Z7]]. In this approach, the original multiwell energy density is replaced by its quasiconvex 
envelop, i. e., by the lowest energy density achievable by the material through the development of 
microstructure. Thus, the determination of the relaxed energy density requires the evaluation of 
all possible microstructures compatible with a prescribed macroscopic deformation. The resulting 
relaxed energy density is quasiconvex p7|], and its minimizers, which represent the macroscopic 
deformation fields, are devoid of microstructure and, thus, more readily accessible to numerical 
methods. In essence, the use of relaxed energy densities in macroscopic boundary value problems 
constitutes a multiscale approach in which the development of microstructure occurs — and is dealt 
with — at the subgrid level. The central problem in this approach is to devise effective means of 
determining the relaxed energy density and of integrating it into macroscopic calculations. 

Unfortunately, no general algorithm for the determination of the quasiconvex envelop of an 
arbitrary energy density is known at present. A fallback strategy consists of the consideration 
of special microstructures only, inevitably resulting in a partial relaxation of the energy density. 
For instance, attention may be restricted to microstructures in the form of sequential laminates 
[ [31] , |36| , [391 , TO [25[ , [24| ]. The lowest energy density achievable by the material through sequential 
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lamination is known as the rank-one convexification of the energy density. Many of the microstruc- 
tures observed in shape-memory alloys [ j22| ] and in ductile single crystals [ j39[ ] may be interpreted 
as instances of sequential lamination, which suggests that the rank-one convexification of the en- 
ergy coincides — or closely approximates — the relaxed energy for these materials. Ductile single 
crystals furnish a notable example in which the rank-one and the quasiconvex envelops are known 
to coincide exactly [Q] . 

In this paper we present a practical algorithm for partially relaxing multiwell energy densities. 
The algorithm is based on sequential lamination and, hence, at best it returns the rank-one con- 
vexification of the energy density. Sequential lamination constructions have been extensively used 
in both analysis and computation [|31], H [I3J, [2^, [39|, [34j |3|, ffl. All microstructures generated 
by the algorithm are in static and configurational equilibrium. Thus, we optimize all the interface 
orientations and variant volume fractions, with the result that all configurational forces and torques 
are in equilibrium. We additionally allow the variants to be arbitrarily stressed and enforce traction 
equilibrium across all interfaces. 

The proposed lamination construction is constrained in an important respect: during a defor- 
mation process, we require that every new micro structure be reachable from the preceding mi- 
crostructure along an admissible transition path. The mechanisms by which microstructures are 
allowed to effect topological transitions are: branching, i. e., the splitting of a variant into a rank- 
one laminate; and pruning, consisting of the elimination of variants whose volume fraction reduces 
to zero. Branching transitions are accepted provided that they reduce the total energy, without con- 
sideration of energy barriers. By repeated branching and pruning microstructures are allowed to 
evolve along a deformation process. The continuation character of the algorithm furnishes a simple 
model of metastability and hysteresis. Thus, successive microstructures are required to be 'close' 
to each other, which restrict the range of microstructures accessible to the material at any given 
time. In general, this restriction causes the microstructures to be path-dependent and metastable, 
and the computed macroscopic response may exhibit hysteresis. 

The proposed relaxation algorithm may effectively be integrated into macroscopic finite-element 
calculations at the subgrid level. We demonstrate the performance and versatility of the algorithm 
by means of a numerical example concerned with the indentation of an Cu-Al-Ni shape memory 
alloy [ |22| ] by a spherical indenter. The calculations illustrate the ability of the algorithm to generate 
complex microstructures while simultaneously delivering the macroscopic response of the mate- 
rial. In particular, the algorithm results in force-depth of indentation curves considerably softer 
than otherwise obtained by direct energy minimization. 

2 Problem Formulation 

Let VL E R 3 be a bounded domain representing the reference configuration of the material. Let 
y(x) : £1 — > R 3 be the deformation and F(x) = Dy(x) be the corresponding deformation gradient. 
We denote the elastic energy density at deformation gradient F E R 3x3 by W(F). We require 
W{F) to be material frame indifferent, i. e., to be such that W(RF) = W(F), Vi? E SO (3) 
and F E R 3x3 . In addition, the case of primary interest here concerns materials such that W(F) 
is not quasiconvex. As a simple example, we may suppose that W(F) has the following special 
structure: Let Wi(F), i = 1, . . . , M be quasiconvex energy densities (see, e. g., [ p3| ] for a definition 
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and discussion of quasiconvexity), representing the energy wells of the material. Then 

W(F) = min W m (F) (1) 

m=0,...,M 

i. e., W(F) is the lower envelop of the functions W m (F). 

A common model of microstructure development in this class of materials presumes that the 
microstructures of interest correspond to low-energy configurations of the material, and that, con- 
sequently, their essential structure may be ascertained by investigating the absolute minimizers of 
the energy. However, the energy functionals resulting from multiwell energy densities such as ([]]) 
lack weak- sequential lower semicontinuity and their infimum is not attained in general [^3]]. The 
standard remedy is to introduce the quasiconvex envelop 

QW(F) = inf / W{F + Du)dx (2) 

\Q\ u6W 1,oo (n) jq 

of W(F), or relaxed energy density. In this expression, W '°°(fl) denotes the space of functions 
whose distributional derivatives are essentially bounded and which vanish on the boundary, and Q 
is an arbitrary domain. Physically, QW(F) represents the lowest energy density achievable by the 
material through the development of microstructure. The macroscopic deformations of the solid 
are then identified with the solutions of the relaxed problem: 

inf { / [QW(Dy) - f ■ y]dx - ! t ■ y\ (3) 

where X denotes some suitable solution space, / is a body force field, t represents a distribution 
of tractions over the traction boundary dil 2 , and the deformation of the body is constrained by 
displacement boundary conditions of the form: 

y = y on dQ x = dtt - dVL 2 (4) 

Thus, in this approach the effect of microstructure is built into the relaxed energy QW(F). The 
relaxed problem defined by QW(F) then determines the macroscopic deformation. 

In executing this program the essential difficulty resides in the determination of the relaxed 
energy QW(F). As mentioned in the introduction, no general algorithm for the determination of 
the quasiconvex envelop of an arbitrary energy density is known at present. A fallback strategy is 
to effect a partial relaxation of the energy density by recourse to sequential lamination [ |3T| , |36| ], 
and the use of the resulting rank-one convexification RW(F) of W(F) in lieu of QW(F) in 
the macroscopic variational problem @. We recall that the rank-one convexification RW(F) of 
W(F) follows as the limit 

RW(F) = lim R k W(F) (5) 

where R W(F) = W(F) and R k W(F) is defined recursively as 

R k W(F) = inf {(1 - X)R k _ 1 W(F - Xa ® N) + XR k ^ x W(F + (1 - A) a ® N), 

X,a,N (6) 

A G [0, 1], a,Ne M 3 , \N\ = l} k > 1 
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Figure 1 : Example of a rank-2 laminate. Ai and A 2 are the volume fractions corresponding to levels 
1 and 2, respectively, and N\ and N 2 are the corresponding unit normals. 

In these expressions, A and 1 — A represent the volume fractions of the /c-level variants, N is the 
unit normal to the planar interface between the variants, and a is a vector (see Fig. []]). 

Unfortunately, a practical algorithm for the evaluation of the rank-one convexification of gen- 
eral energy densities 'on the fly' does not appear to be available at present. Dolzmann [ |25| ] has 
advanced a method for the computation of the exact rank-one convexification of an arbitrary energy 
density in two dimensions. However, extensions of the method to three dimensions do not appear 
to be available at present. In addition, the method requires the a priori tabulation of RW(F) over 
all of M 2x2 , which limits its applicability to large-scale computing. 

An additional complication arises from the fact that the cyclic behavior of martensitic materials 
often exhibits hysteresis. Under these conditions, the response of the material is path-dependent 
and dissipative, and, therefore, absolute energy minimization does not furnish a complete model 
of material behavior. The modeling of hysteresis requires consideration of entire deformation 
processes, rather than isolated states of deformation of the material. A framework for the un- 
derstanding of hysteresis may be constructed by assuming that the evolution of micro structure is 
subject to a continuity requirement, namely, to the requirement that successive microstructures be 
close to each other in some appropriate sense. This constraint restricts the range of microstructures 
which the material may adopt at any given time and thus results in metastable configurations. The 
particular sequence of metastable configurations adopted by the material may be path dependent, 
resulting in hysteresis. The connection between metastability and hysteresis has been discussed by 
Ball, Chu and James [|9p. 



3 A sequential lamination algorithm 

The problem which we address in the remainder of this paper concerns the formulation of ef- 
ficient algorithms for the evaluation of RW(F), and extensions thereof accounting for kinetics 
and nonlocal effects, with specific focus on algorithms which can be effectively integrated into 
large-scale macroscopic simulations. We begin by reviewing basic properties of sequential lami- 
nates for subsequent reference. More general treatments of sequential lamination may be found in 

[0[^Q[i^g§[3ejp. 

Uniform deformations may conventionally be categorized as rank- zero laminates. A rank-one 
laminate is a layered mixture of two deformation gradients, F~ , F + G IR 3x3 . Compatibility of 
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deformations then requires F^ to be rank-one connected, i. e., 

F + -F =a®N (7) 

where oel 3 , and N £ M 3 , \N\ = 1, is the normal to the interface between the two variants of 
deformation. Let A 1 * 1 , 

A- + A+ = l, A ± e[0,l], (8) 

denote the volume fractions of the variants. Then, the average or macroscopic deformation follows 
as 

F = \~~F~ + X + F + (9) 
If F and {a, X ± , N} are known, then the deformation in the variants is given by: 

F + = F + \-a®N 

F =F - \ + a®N. 1 j 

and, thus, F and {a, X ± , N} define a complete set of-deformation and configurational-degrees of 
freedom for the laminate. Following Kohn pl[], a laminate of rank- A; is a layered mixture of two 
rank-(/c — 1) laminates, which affords an inductive definition of laminates of any rank. As noted 
by Kohn [fTfl], the construction of sequential laminates assumes a separation of scales: the length 
scale l k of the A;th-rank layering satisfies l k <C lk-i- 

Sequential laminates have a binary-tree structure. Indeed, with every sequential laminate we 
may associate a graph G such that: the nodes of G consist of all the sub-laminates of rank less 
or equal to the rank k of the laminate; and joining each sub-laminate of order 1 < I < k with 
its two constituent sub-laminates of order / — 1. The root of the graph is the entire laminate. 
Two sequential laminates will be said to have the same structure (alternatively, topology or layout) 
if their graphs are identical. Evidently, having the same structure defines an equivalence relation 
between sequential laminates, and the set of all equivalence classes is in one-to-one correspondence 
with the set B of binary trees. 

Let i = 1, . . . , n be an enumeration of the nodes of G. Then, to each node i we may associate 
a deformation Fi. The root deformation is the average or macroscopic deformation F. Each node 
in the tree has either two children or none at all. Nodes with a common parent are called siblings. 
Nodes without children are called leaves. Nodes which are not leaves are said to be interior. The 
deformations of the children of node i will be denoted Ff. Each generation of nodes is called a 
level. The root occupies level of the tree. The number of levels is the rank k of the tree. Level 
I contains at most 2 l nodes. The example in Fig. ^ represents a rank-two laminate of order four. 
The three leaves of the tree are nodes F + , F~ + and F . The interior nodes are F and F~ . The 
children of, e. g., node F are nodes F~ + and F . 

Compatibility demands that each pair of siblings be rank-one connected, i. e., 

F+ -FT = a t ®N u iel G (11) 

where 6 M 3 , Ni e M 3 , \N\ = 1, is the normal to the interface between F~ and Ff, and X G 
denotes the set of all interior nodes. Let \f, 



Ar + A+ = 1, A±G[0,1], 



(12) 
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Figure 2: Example of a rank-2 laminate. In this example, V\ = A A + and U2 = X X 



denote the volume fractions the variants Ff within node i. Then, the deformation of the parent 
variant is recovered in the form 

Ft = X-FT + A+F+ (13) 

If Fi and {a^ \ f, iVj} are known for an interior node i, then the deformation of its children is 
given by: 

Ff = Fi+ Xra^Ni 
FT = Fi — At a, ® Nt. 

Therefore, F and {aj, Xf, iVj, i El} define a complete and independent set of degrees of freedom 
for the laminate. A recursive algorithm for computing all the variant deformations Fi, i = 1, . . . , n 
from F and {ctj, Xf, N i: , i gIg} has been given by pQ[]. 

We shall also need the global volume fractions z// of all leaves I E Cq, where Cq denotes the 
collection of all leaves of G. These volume fractions are obtained recursively from the relations: 



vf = Xfv h i E 1 G (15) 
with z/root = 1 for the entire laminate, and satisfy the relation: 

J2 ui = 1 (16) 

Thus, v\ represents the volume occupied by leaf I as a fraction of the entire laminate. The Young 
measure (e. g., [|3~7|]) of the laminate consist of a convex combination of atoms 5 Ft (F) with weights 
v u I E C. 

The average or macroscopic stress of the laminate may be expressed in the form 



(17) 
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where 

P, = W, F (F,), leC G (18) 

are the first Piola-Kirchhoff stresses in the leaves. The average stress may be computed by recur- 
sively applying the averaging relation: 

Pi = \~Pr + X+P+, % e 1 (19) 

starting from the leaves of the laminate. A recursive algorithm for computing the average stress P 
from {Pi, I G £} and {A^, % G Z G } has been given by [|40|]. 

3.1 Microstructural equilibrium 

We begin by investigating the mechanical and configurational equilibrium of sequential laminates 
of given structure. Thus, we consider an elastic material of strain-energy density W(F) subject to 
a prescribed macroscopic deformation F G M 3x3 . In addition, we consider all sequential laminates 
with given graph G. 

The equilibrium configurations of the laminate then follow as the solutions of the constrained 
minimization problem: 



GW(F) = inf > ViW{F{) (20) 

Afe[0,l], ieI G (21) 
|iVi| = l, tela (22) 



where the Fi are obtained from the recursive relations dl4]). The effective or macroscopic energy of 
the laminate is GW(F). If W(F) is quasiconvex then it necessarily follows that a, = 0, Vz G T G , 
and GW(F) = W(F). 

It is interesting to verify that the solutions of the minimization problem d20l ) are in both force 
and configurational equilibrium. Thus, assuming sufficient smoothness, the stationarity of the 
energy with respect to all deformation jump amplitudes yields the traction equilibrium equations: 

(Pt-PT)-Ni = o, iel G (23) 

Stationarity with respect to all normal vectors yields the configurational-torque equilibrium equa- 
tions: 

[a t ■ (Pt - Pr)] xNi = 0, i G Tq (24) 

Finally, stationarity with respect to all volume fractions yields the configurational-force equilib- 
rium equations: 

f l = (W t + -Wr)-(KPt + KP7)-(a i ®N i ) = 0, i G Tq (25) 



where fi is the configurational force which drives interfacial motion [|14|]. It bears emphasis that 
the leaf deformations Fi may in general be arbitrarily away from the minima of W(F), and thus 
the equilibrium equations ( p3| ) must be carefully enforced. In addition, the minimization ( pQj ) 
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has the effect of optimizing the volume fractions of all variants and the corresponding interface 
orientations. If all the preceding stationarity conditions are satisfied, then it is readily verified that 



the average or macroscopic deformation (17) is recovered as 



p = GW, F (F) (26) 

which shows that GW(F) indeed supplies a potential for the average or macroscopic stress of the 
laminate. 

The case in which the energy density W(F) possesses the multiwell structure (|T]) merits special 
mention. In this case, the minimization problem (|20l - [22l ) may be written in the form: 



GW(F)= inf X^W^F,) (27) 

\ ,N U i e 1 G } i eCo 
{ mi e {l,...,M}, l g C G } 

Xf g [0, 1], iel G (28) 
|JVi| = i, «el G (29) 

where m = {mi G {1, . . . , M}, I G C G } denotes the collection of wells which are active in each 
of the leaves. This problem may be conveniently decomposed into two steps: a first step involving 
energy minimization for a prescribed distribution of active wells, namely, 

G m W(F)= inf \> z W m! (F z ) (30) 

A±G[0,1], %EX G (31) 
1^1 = 1, iel G (32) 

followed by the optimization of the active wells, i. e., 

GW(F) = inf G m W(F) (33) 

{ mi e{i,...,M}, i&C G } 



It should be carefully noted that the minimizers of problem ( |20| - [22| ) may be such that one or 
more of the volume fractions Xf take the limiting values of or 1. We shall say that a graph G 
is stable with respect to a macroscopic deformation F if at least one minimizer of (^(J-^2|) is such 
that 

A±G(0,1), Vi G Tq (34) 

and we shall say that the graph is unstable or critical otherwise. The presence of sub-trees of 
zero volume in an unstable graph is an indication that the graph is not 'right' for the macroscopic 
deformation F, i. e., the graph is unable to support a nontrivial micro structure consistent with F. 
Unstable graphs are mathematically contrived and physically inadmissible, and, as such, should be 
ruled out by some appropriate means. This exclusion may be accomplished, e. g., by the simple 
device of assigning the offending solutions an infinite energy, which effectively rules them out 
from consideration; or by defining solutions modulo null sub-trees, i. e., sub-trees of vanishing 
volume. In the present approach, we choose to integrate the exclusion of null sub-trees into the 
dynamics by which microstructures are evolved, as discussed next. 
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3.2 Microstructural evolution 

The problem ( |20l - [22l ) may be regarded as a partial rank-one convexification of W(F) obtained by 
prescribing the graph G of the test laminates. The full rank-one convexification follows from the 
consideration of all possible graphs, i. e., 

RW(F) = inf GW(F) (35) 

where, as before, B is the set of all binary trees. In the particular case of energy densities of the 
form ((!]), we alternatively have 

RW(F) = inf G m W(F) (36) 

G e B 

{ mi G {1, . . . , M}, I G C G } 
It is clear that this problem exhibits combinatorial complexity as the rank of the test laminates 



increases, which makes a direct evaluation of (35) or (B6J) infeasible in general. 



Problems of combinatorial complexity arise in other areas of mathematical physics, such as 
structural optimization and statistical mechanics. Common approaches to the solution of these 
problems are to restrict the search to the most 'important' states within phase space, or importance 
sampling; or to restrict access to phase space by the introduction of some form of dynamics. In 
this latter approach the states at which the system is evaluated form a sequence, or 'chain', and 
the next state to be considered is determined from the previous states in the chain. If, for instance, 
only the previous state is involved in the selection of the new state, a Markov chain is obtained. In 
problems of energy minimization, a common strategy is to randomly 'flip' the system and accept 
the flip with probability one if the energy is reduced, and with a small probability if the energy is 
increased. 

In other cases, the system possesses some natural dynamics which may be exploited for compu- 
tational purposes. A natural dynamics for problem ( Tj"5j ) may be introduced as follows. Evidently, 
the relevant phase space for this problem is B, the set of all binary trees, and the aim is to define 
a flow G(t) in this phase space describing the evolution of the micro structure along a deformation 
processes F(t). Here and subsequently, the real variable t > denotes time. A natural dynamics 
for G(t) is set by the following conditions: 

1. G(t) must be stable with respect to F(t). 

2. G(t) must be accessible from G(t~) through a physically admissible transition. 

The first condition excludes laminates containing null subtrees, i. e., subtrees of zero volume. The 
second criterion may be regarded as a set of rules for microstructural refinement and unrefinement. 
In order to render these criteria in more concrete terms, we adopt an incremental viewpoint 

and seek to sample the micro structure at discrete times to — 0, . . . , t n , t n+ i, Suppose that 

the micro structure is known at time t n and we are given a new macroscopic deformation F n+1 = 
F(t n+ i). In particular, let G n be the graph of the micro structure at time t n . We consider two 
classes of admissible transitions by which a new structure G n +i may be reached from G n : 



1. The elimination of null subtrees from the graph of the laminate, or pruning. 
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2. The splitting of leaves, or branching. 



Specifically, we refer to as branching the process by which a leaf is replaced by a simple laminate. 
The criterion that we adopt for accepting or rejecting a branching event is simple energy minimiza- 
tion. Thus, let I G £ Gn be a leaf in the microstructure at time t n , and let F™ be the corresponding 
deformation. The energetic 'driving force' for branching of the leaf I may be identified with 

ff = W(F?) - RiW(F™) (37) 

where RiW(F) is given by We simply accept or reject the branching of the leaf I G C Gn 
according as to whether /" > or /" < 0, respectively. 

In the particular case in which W(F) is of the form (P, the evaluation of R 1 W(F) may be 
effected by considering all pairs of well energy densities, one for each variant. However, since 
the well energy densities W m (F) are assumed to be quasiconvex and we rule out variants of zero 
volume, we may exclude from consideration the cases in which both variants of the laminate are in 
the same well. The evaluation of R\W{F) is thus reduced to the consideration of all distinct pairs 
of well energy densities. 

The precise sequence of steps followed in calculations are as follows: 



1. Initialization: Input F n+1 , set G n+ \ = G n . 

2. Equilibrium: Equilibrate laminate with G n +i held fixed. 

3. Evolution: 



(a) Are there null subtrees? 

i. YES : Prune all null subtrees, GOTO (2). 

ii. NO : Continue. 

(b) Compute all driving forces for branching {f? +1 , I G £ Gn+1 }. Let /;™tx = max 'e£ Gn 1 fi +l - 
Is /, n+1 > 0? 

** 'max 

i. YES : Branch leaf Z max , GOTO (2). 

ii. NO : EXIT. 



Several remarks are in order. The procedure just described may be regarded as a process of 
continuation, where the new microstructure is required to be close to the existing one in the sense 
just described. Evidently, since we restrict the class of microstructures which may arise at the end 
of each time step, there is no guarantee that this continuation procedure deliver the solution of ( |35| ) 
for all t > 0. However, metastability plays an important role in many systems of interest, and the 
failure to deliver the absolute rank-one convexification at all times is not of grave concern in these 
cases. Indeed, the continuation procedure described above may be regarded as a simple model of 
metastability. 

In this regard, several improvements of the model immediately suggest themselves. Thus, 
the branching criterion employed in the foregoing simply rules out branching in the presence of 
an intervening energy barrier, no matter how small, separating the initial and final states. An 
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improvement over this model would be to allow, with some probability, for transitions requiring 
an energy barrier to be overcome, e. g., in the spirit of transition state theory and kinetic Monte 
Carlo methods. However, the implementation of this approach would require a careful and detailed 
identification of all the paths by which branching may take place, a development which appears 
not to have been undertaken to date. 



4 Illustrative examples 

As a first illustration of the sequential lamination algorithm presented in the foregoing, we apply 
it to a simple model of a Cu-Al-Ni shape-memory alloy, a material which has been extensively 
investigated in the literature (cf [|7], [8|, [22|, ^6[, and references therein). Photomicrographs taken 
from the experiments of Chu and James [ |2"2"| ] (also reported by P^]) reveal sharp laminated mi- 
crostructures, often of rank two or higher. In order to exercise the algorithm, in the examples that 
follow we simply take the material through a prescribed macroscopic deformation path. 



4.1 Material model 

Cu-Al-Ni undergoes a cubic to orthorhombic martensitic transformation at around room temper- 
ature and has, therefore, six variants in the martensitic phase. The deformation undergone by the 
material in transforming from austenite to an unstressed variant of martensite may be described by 
a stretch tensor U m , m = 1, . . . 6. For Cu-Al-Ni , these are [E2L 5TJ]: 
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(38) 



(39) 



where £ = 1.0425, rj = 0.0194 and ( = 0.9178, and all components are referred to the cubic 
axes of the austenitic phase. For purposes of illustration of the sequential lamination algorithm we 
adopt a simple energy density of the form ([]]), with well energy densities 



W n (F) 



F'F-I) C (F 1 F- I) 



for the austenitic phase, and 
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(41) 
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for the martensitic phases. In these expressions C r 
bottom of the variants. These are (in MPa): 
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(43) 



where Cn = 189, C 22 = 141 and C33 = 205, C 12 = 124, C 13 = 45.5, C 23 = 115, C 44 = 54.9, 
C55 = 19.7, 6*66 = 62.6 and with the moduli of the remaining martensitic variants following 
by symmetry. The energy density defined by ([]]) and ( ]4T| ) is material frame indifferent, results 
in stress-free states at the bottoms of all wells, i. e., at F = I and F = RU m , R E SO(3), 
m = 1, . . . , 6, assigns equal energy density to all unstressed variants, and exhibits all the requisite 
material symmetries. 



4.2 Optimization 

In the examples presented here, and in the finite element calculations presented in the follow- 
ing section, problem ( |2"0] - |2~2"| ) is solved using Spellucci's sequential quadratic programming (SQP) 
algorithm for constrained minimization. The SQP algorithm is an interative procedure which re- 
quires an initial guess in order to start the iteration. In calculations we begin by setting the initial 
values of {a*, \ iVj, i e Xc n+1 } at time t n+ i equal to the converged values at time t n . The main 
issue arises in evaluating possible branching events, as in this case new interfaces arise for which 
no previous geometrical information exists. The selection of initial guesses for <2j and A^, i E X G , 
offers no difficulty, and we simply set a« = and = or 1 . The choice of the initial value of the 
new interface normals Ni requires more care since it strongly biases the resulting micro structure. 

We have investigated two ways of initializing the normals Ni arising during branching. A first 
approach is based on sampling the unit sphere uniformly. Thus, we simply select initial values 
of Ni uniformly distributed over the unit sphere with some prespecified density and select the 
solution which results in the least energy. If two or more branched configurations possess the same 
energy, we select one at random. This exhaustive search approach is effective but costly owing to 
the large number of cases which need to be considered. 

The second approach consists of priming the iteration using an initial guess derived from Ball 
and James [0, [8]] constrained theory for shape-memory materials. In this theory, the elastic moduli 
are presumed large compared to the transformation stresses, so that the geometry of the laminate 
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Table 1: Vector e arising in the twinning relations for Cu-Al-Ni JT7|] . The vectors {ei,e 2 ,e 3 } 
correspond to the cube directions in the austenite phase. 



can be obtained, to a first approximation, directly from the transformation strains. Conveniently, 
all resulting twinning relations between every distinct pair of martensitic wells can be tabulated 
beforehand. For Cu-Al-Ni this tabulation has been carried out by Bhattacharya, Li and Luskin [|P7j]. 
As expected from general theory, the only non-trivial rank-one connections take place between 
distinct variants of martensite. Specifically, we seek Q G £0(3), a, N G IR 3 , |JV| = 1, such that 

QU m = U n + a<g)N, m,n = l,...,6,m^n (44) 



There are two solutions to this equation. Let e = e/\e\, with the vector e as in Table |]. The first 
solution is 



a 



and the second solution is 
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(45) 
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1, . . . , 6, n 7^ m (46) 



These results from the constrained theory may conveniently be used to start a branching calculation 
when the energy density is of the form ([]]). In this case, the branching calculation entails the 
consideration of every pair of well energy densities W m (F) and W n (F), m, n = 1, . . . , 6, m ^ n, 
in the new leaves. The attendant iteration may then be started from the constrained solutions 
{a, N} described above, with the initial value of A ± determined by means of a line search. 



4.3 Martensite-martensite transition 

Our first test case concerns the macroscopic deformation process 

F(t) = (1 - t)U x + tU 2 , t G [0, 1] (47) 

which takes the material from one variant of martensite to another. Fig. [5] shows the evolution of 
the energy, the component P i3 of the first Piola-Kirchhoff stress, and the volume fraction A of U 2 , 
respectively, for the unrelaxed and relaxed cases. In this calculations, the branching constructions 
employ the constrained geometry as initial guess, as discussed in the foregoing. 

The unrelaxed response shown in Fig. ^exhibits an abrupt transition from the initial to the final 
variant, as no mixed states are allowed to develop during the deformation process. By contrast, the 
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Figure 3: Martensite-to-martensite transition example. Comparison of unrelaxed (left) and relaxed 
(right) energies, stresses, and volume fractions. 
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step step 



(a) (b) 

Figure 4: Simple shear example, a) Pi 2 component of the first Piola-Kirchhoff stress tensor as 
a function of applied shear strain. The graph of the microstructures predicted by the relaxation 
algorithm are also shown inlaid, b) Unrelaxed and relaxed energies during loading. 



relaxation algorithm results in the development of a rank-one laminate immediately following the 
onset of deformation. It should be carefully noted that the material is allowed to develop laminates 
of arbitrary rank, and that the persistency of a rank-one laminate is due to the fact that both leaves 
are stable against branching. The computed volume fraction A increases linearly from to 1. As a 
result of this evolving microstructure, the energy of the material is fully relaxed, and the material 
remains unstressed. In this example, the unloading response exactly traces in reverse the loading 
response, and hence no hysteresis is recorded. 

4.4 Simple shear 

Our second test case concerns macroscopic simple shear on the plane (010) and in the direction 
(100). The material is taken to be initially undeformed, and the shear deformation is initially 
increased from zero up to a maximum value and then decreased to back zero. The calculations are 
carried out excluding the austenitic well from the definition ([]]) of the energy and considering the 
six martensitic wells only. 

Fig. 0b shows a comparison of the computed unrelaxed and relaxed energies. By the exclusion 
of the austenitic well the material is forced to develop a rank-five laminate in its initial undeformed 
configuration. The graph of this laminate, and of all laminates which subsequently arise during 
deformation, is shown inlaid in Fig. |]a. As is evident from Fig. |]b, the relaxation algorithm 
ostensibly succeeds at fully relaxing the energy of the material. With increasing deformation, the 
computed microstructure undergoes transitions to rank-four and three laminates. A first rank-three 
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laminate of order thirteen is first predicted which subsequently simplifies to a rank-three laminate 
of order seven. As a result of this micro structural evolution, the relaxed energy remains well below 
the unrelaxed energy through the deformation, Fig. |]b. Upon unloading, the order-seven rank- 
three laminate is maintained down to zero deformation, suggesting that the initial microstructures 
computed during loading are metastable. Indeed, the unloading stress-strain curve lies below the 
loading one, resulting in a certain amount of hysteresis, which suggests that the order-seven rank- 
three laminate is indeed more efficient that the precursor microstructures. 



5 Nonlocal extension 



The simple branching criterion (p7[), which accounts for the energies of the variants only, neglects 
any energy barriers as may opposed the transformation and may lead to runaway refinement of 
the micro structure. In order to limit branching and, additionally, to estimate the size of the mi- 
crostructure, we follow Ball and James [Q] and take into consideration two additional sources of 
energy: the energy of the twin boundaries, and the mismatch energy contained within the boundary 
layers separating pairs of leaves. For instance, Boullay et al. [JTSp, and James et al. pD[ ] have in- 
vestigated in detail the structure of branched needle microstructures that develop within the misfit 
boundary layers, e. g., at the edge of a martensite laminate. Such level of detail is well beyond the 
scope of this work. Our aim here is to derive a rough estimate of the misfit energy amenable to a 
straightforward calculation. 



Sampling Points 




Boundary Layer 



Figure 5: Schematic of interpolation boundary layer and scheme used to estimate misfit energy. 



One such simple estimate may be derived as follows. Begin by enforcing 'rigid-device' bound- 
ary conditions 

y BL = y + F{x-x ), x 3 = (48) 

where F is the average deformation in the laminate, x is a material point within the reference 
configuration of the laminate, and y = y(x ) is its position on the deformed configuration. This 
insulates the laminate from the details of the adjacent deformation field in the regions x 3 > and 
£3 < — I. A simple interpolating deformation mapping is then: 

y BL (aO = y - (y(x) - y )^ + F(x - x ) (l + ^ ) , -A < x 3 < (49) 
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where y(x) is the deformation mapping of the laminate. The boundary layer at x 3 = —I can be 
given an identical treatment. The corresponding deformation mapping is 



-iBL, 



x) = -F(aj)^ + F (l + ^) + ±[F(x - x ) - (y(x) - y )} ® e 3 (50) 



We proceed to estimate the elastic energy of the region defined by the intersection of each variant 
with the boundary layer by a simple one-point quadrature rule. Let x ± be a pair of consecutive 
sampling points, Fig. |5], chosen such that xf = —A/2, and select x = x for simplicity. Then, it 
follows immediately from (^0|) that 

f bl- = F BL (a:-) = -(F~ + F) (51) 



Likewise, 



F BL+ = F Bh (x+) = + F) + [F(x + - X-) - (y + - y")] ® e 3 (52) 



But, since y(a;) is piecewise linear, it follows that 

y+ -y = F-(\-(x + - x-)) + F + (A + (> + - aT)) = F(a; + - a; - ) (53) 

and 

F BL + = I(F+ + F) (54) 
2 

Collecting the above results we finally have 

F BL± = ^(F ± + F) (55) 

Within this approximation, the misift boundary-layer energy density finally evaluates to 

W BL = X~[W(F Bh ~) - W(F-)} + X + [W(F BL+ ) - W(F + )} (56) 
which furnishes a remarkably simple (though rough) estimate. We note that 

F BL+ - F BL - = ^{F + - F~) = ^a® N (57) 

and 

\-F BL - + A + F BL+ = F (58) 

Thus, the deformations F BL± axe rank-one compatible and match the average deformation of the 
laminate. Since the deformations F minimize the energy of the laminate among all rank-one 
laminates with average deformation F, it follows that W BL > 0, i. e., W BL does indeed represents 
an excess energy. 

For simplicity, we assume that the twin-boundary energy T per unit area is a constant indepen- 
dent of the deformation of the variants. Combining the preceding estimates, it follows that total 
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excess or nonlocal energy due to the twin boundaries and the misfit boundary layers contained 
within a region of the laminate of dimensions L x L x I is: 



T 2A 

F + T 

Taking A = / c /2, for definiteness, this expression reduces to 



£NL = L * l )h- + ^±W Bh } (59) 



£NL = hH I t_ + ' W BL . (6Q) 



r i c 
F + T 

This excess energy may now be minimized with respect to l c , with the result: 



' c = Vi^ <61) 

which affords an estimate of l c . The corresponding minimum excess energy per unit volume is 



W " L = In = 2 vT~ <«> 

We note that this excess energy grows as Z -1 / 2 , which tends to suppress microstructural refinement. 
In calculations, we interpret the excess energy density W NL as an energy barrier for branching. 
Consideration of this energy barrier has the effect of reducing the local branching driving force 
(0)to 

ft = W(F?) - R 1 W(F] L ) - W NL (F") (63) 

which effectively introduces a lower cut-off for the size of the microstructure and eliminates the 
possibility of runaway microstructural refinement. 



6 Application to the finite-element simulation of indentation in 
Cu-Al-Ni 

The sequential lamination algorithm developed in the foregoing may conveniently be taken as a 
basis for multiscale simulation in situations in which there is a strict separation of scales: a macro- 
scopic scale characterized by slowly-varying smooth fields; and a much smaller scale commensu- 
rate with the size of the evolving microstructure. As remarked by several authors [ PD[ , ^6] , |T9"| , [23[ ], 



these problems may be solved effectively by pushing the microstructure to the sub-grid scale, while 
solving the well-posed relaxed problem on the computational grid. 

In this section we present selected examples of application of this multiscale approach in which 
the macroscopic problem is solved by the finite-element method, while the effective behavior is 
computed, simultaneously with the macroscopic solution, at the Gauss-point level using the se- 
quential lamination algorithm developed in the foregoing. The particular problem considered con- 
cerns the quasistatic normal indentation of a Cu-Al-Ni shape-memory alloy by a spherical indenter. 
The domain of analysis and the computational mesh are shown in Fig. ||. The analysis is reduced to 
one quarter of the entire domain for simplicity. In particular, solutions exhibiting broken symmetry 
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are ruled out by the analysis. The size of the computational domain is 20mm x 20mm x 20mm. 
The radius of the indenter is 15mm. The specimen is fully supported over its entire base, and the 
remainder of its boundary is free of tractions. The computational mesh contains 254 nodes and 
105 ten-node quadratic tetrahedral elements. Contact between the indenter and the specimen is as- 
sumed to be frictionless and is enforced by a penalty method To ensure that the jacobian J of 



the deformation remains positive in all variants at all times, the energy of each well is augmented 
by a term of the form [PSR 



w-\j) = r^ ^ 'it (64) 



C(J 2 + J- 2 - 2) 2 , J < 1 
0, otherwise 

where C is a constant chosen sufficiently small to minimize the effect on the total energy. By 
design, W vol (J) and its first and second derivatives vanish at J = 1. In addition, the twin-boundary 
and misfit energies are accounted for as part of the branching criterion as a means of introducing 
a lower cutoff for the laminate size and preventing runaway microstructural refinement. In all 
calculations, the twin-boundary energy per unit area T is set to 1 J/m 2 . The maximum size of the 
laminate at a particular Gauss point is set to the element size, in keeping with the assumption that 
the laminate accounts for sub-grid structure in the solution only, and that coarser structures are 
adequately resolved by the mesh. 

The finite-element solution is obtained by dynamic relaxation followed by a preconditioned 
conjugate-gradient iteration [03]]. The high level of concurrency in the constitutive calculations 



was exploited via an MPI-based parallel implementation Q42Q on the ASCI Blue multiprocessing 



computer. Performance studies showed excellent load balancing and scalability. 



Fig. [7] shows the unrelaxed deformed configurations, and the corresponding distribution of 
active energy wells at the Gauss points of the mesh, at depths of indentation of 0. 150 and 0.375 mm. 
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Figure 7: Cross-sections and energy-density contours for two unrelaxed solutions at depths of 
indentation: a) 0.150 mm, and b) 0.375 mm. The symbols designate the energy well which is 
activated at each Gauss point of the mesh. 
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As is evident from this figure, two energy wells become active during indentation. The transformed 
zone under the indenter grows with depth of indentation, but the fineness of the variant arrangement 
is severely limited by the mesh size. Correspondingly, the total energies and indentation forces 
recorded during indentation are comparatively high, Fig. |9|. In this figure the energy has been 
normalized by E = y C^ stemte , where Vq is the volume of the undeformed specimen, while the 
force has been normalized by F = E / R Indenter . 

The relaxed solution obtained using the sequential lamination algorithm differs markedly from 
the unrelaxed solution just described, Fig. |8|. Thus, the relaxed deformation field is accompanied 
by the development of well-defined microstructures at the subgrid level. Some of the laminates 
generated by the sequential lamination algorithm are quite complex, reaching rank two. Of note 
is the appearance of a de-twinned zone directly under the indenter. The effect of relaxation on 
the total energy and indentation force is quite marked, Fig. |9|, with the relaxed values lying well 
below the unrelaxed ones. Unloading exhibits the path-dependent nature of the algorithm, with 
the micro structure established at maximum load remaining in place during much of the unloading 
process, which in turn results in a soft response. The fineness of the microstructure is somewhat 
overpredicted by the calculations, with some of the variants attaining sub-micron thicknesses. In 
view of (|6T|), this excessive refinement may owe to a low value of the twin-boundary energy T, or 
to an overestimation of the misfit energy W BL , or both. 



7 Summary and Concluding Remarks 

We have presented a practical algorithm for partially relaxing multiwell energy densities. The 
algorithm is based on sequential lamination, but it is constrained in such a way that successive 
microstructures occurring along a deformation path are close to each other in a certain sense: the 
new microstructure should be reachable from the preceding one by a combination of branching and 
pruning operations. All microstructures generated by the algorithm are in static and configurational 
equilibrium. In particular, we optimize all the interface orientations and variant volume fractions, 
with the result that all configurational forces and torques are in equilibrium. We additionally allow 
the variants to be arbitrarily stressed and enforce traction equilibrium across all interfaces. Owing 
to the continuity constrained imposed upon the micro structural evolution, the predicted material 
behavior may be path-dependent and exhibit hysteresis. 

In cases in which there is a strict separation of micro and macrostructural lengthscales, the 
proposed relaxation algorithm may effectively be integrated into macroscopic finite-element cal- 
culations at the subgrid level. We have demonstrated this aspect of the algorithm by means of a 
numerical example concerned with the indentation of an Cu-Al-Ni shape memory alloy [ ]22| ] by 
a spherical indenter. The calculations illustrate the ability of the algorithm to generate complex 
microstructures, resulting in force-depth of indentation curves considerably softer than otherwise 
obtained by direct energy minimization. 

Several improvements of the present approach immediately suggest themselves. The relaxation 
algorithm, in its present form, does not account for energy barriers for branching. Thus, a signif- 
icant improvement over this model would be to permit, with probability less than one, transitions 
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Indentation = 0.15 Indentation = 0.375 



(a) (b) 

Figure 8: Cross-section and energy-density contours for relaxed solution at an indentation depth of: 
a) 0.150 mm, and b) 0.375 mm. The symbols indicate the rank of the microstructure at the Gauss 
points. The insets depict the geometry of the microstructure at the indicated sampling points, with 
each color representing an individual well, and are of identical size oriented such that the left face 
corresponds to the cross-section plane. 
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Figure 9: Normalized total energy vs. normalized depth of indentation; b) Normalized indentation 
force vs. normalized depth of indentation. 

requiring an energy barrier to be overcome, e. g., in the spirit of transition state theory and kinetic 
Monte Carlo methods. This extension would require a careful and detailed identification of all 
the paths by which branching may take place, and the attendant energy barriers. Another signif- 
icant improvement would be to relax the configurational force equilibrium constraint and replace 
it by a kinetic relation governing interfacial motion [Jl], ^, [14]]. Kinetic relations of this form can 
effectively be integrated into the variational principle with the aid of time discretization [^9] ftTj ]. 
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